{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive Derivative-Assembled Pseudo-Trotter ansatz Variational Quantum Eigensolver [1] (ADAPT-VQE)\n",
    "\n",
    "The Unitary Coupled Clusters truncated at the level of double excitations (UCCSD)\n",
    "$$ |\\psi^{UCCSD}\\rangle \\equiv e^{T(1)+T(2)}|\\psi^{HF}\\rangle, \\label{eq:1}\\tag{1}$$\n",
    "where\n",
    "$$ T(1)=\\sum_{ia} \\hat \\tau_{ia} = \\sum_{ia} \\tau_{ia} (a^{\\dagger}_{a} a_{i} - a^{\\dagger}_{i} a_{a}),$$\n",
    "$$T(2) = \\sum_{iajb} \\hat \\tau_{iajb} = \\sum_{iajb} \\tau_{iajb} (a^{\\dagger}_{a}  a^{\\dagger}_{b} a_{j} a_{i} - a^{\\dagger}_{i} a^{\\dagger}_{j} a_{b} a_{a}).$$\n",
    "is a natural ansatz for Variational Quantum Eigensolver since it is parametrized by the unitary operators and has a polynomial number of variational parameters. For the practical application on the quantum computer, it is desired to decompose unitary operator in Eq. (\\ref{eq:1}) into a product of unitary operators acting on a few qubits. This can be achieved via the first-order Trotter decomposition. Because of the variational flexibility of the ansatz (parameters are optimized after the Trotterization), one can sufficiently reproduce the UCCSD results performing only single Trotter step\n",
    "$$ |\\psi^{UCCSD}\\rangle \\approx \\prod_{ia}e^{\\hat \\tau_{ia}} \\prod_{iajb}e^{\\hat \\tau_{iajb}} |\\psi^{HF}\\rangle. \\label{eq:2}\\tag{2}$$\n",
    "The ADAPT-VQE algorithm suggests limiting the number of unitary operators in Eq. (\\ref{eq:2}) to lower the depth of the circuit and make the ansatz even more practical for quantum computing. This is achieved by the systematic extension of the ansatz by one excitation operator at a time, starting from the Hartree Fock state. Specifically, having specified the pool of all possible single and double excitation operators\n",
    "$\\{ \\hat A_1, \\hat A_2, \\dots, \\hat A_n \\}$ ($\\hat A_i \\in \\{ a^{\\dagger}_a a_i - a^{\\dagger}_i a_a; a^{\\dagger}_{a}  a^{\\dagger}_{b} a_{j} a_{i} - a^{\\dagger}_{i} a^{\\dagger}_{j} a_{b} a_{a} \\}$), at first one evaluates the gradient $$\\frac{\\partial E}{\\partial \\theta_m} = \n",
    "\\left.\\frac{\\partial}{\\partial \\theta_m} \\langle \\psi^{HF} |e^{-\\theta_1 \\hat A_1} e^{-\\theta_2 \\hat A_2} \\ldots\n",
    " \\hat H \\ldots e^{\\theta_2 \\hat A_2} e^{\\theta_1 \\hat A_1}| \\psi^{HF} \\rangle \\right|_{ \\forall \\theta_i=0} \n",
    " =\\langle \\psi^{HF} | [ \\hat H,\\hat A_m] | \\psi^{HF} \\rangle . \\label{eq:3}\\tag{3}$$\n",
    " Then ansatz is extended with the excitation operator,$\\hat A_k$, corresponding to the highest element of the gradient. Afterwards, all or only the new variational parameters are optimized. Continuing this process until the energy, the L$^2$ norm of the gradient or the maximum element is converged one recovers ADAPT-VQE [1]\n",
    "$$|\\psi^{ADAPT}\\rangle \\equiv (e^{\\hat \\tau_N}) (e^{\\hat \\tau_{N-1}}) \\ldots (e^{\\hat \\tau_2}) (e^{\\hat \\tau_1}) |\\psi^{HF}\\rangle, \\label{eq:4}\\tag{4}$$\n",
    "where $\\hat \\tau_p$ can be either $\\hat \\tau_{ia}$ or $\\hat \\tau_{iajb}$. It should be emphasized that the operator  pool is never drained and always contains all possible single and double excitation operators, thus allowing some excitation operator,$\\hat A_k$, enter the ansatz multiple times.\n",
    "\n",
    "[1] Harper R. Grimsley, Sophia E. Economou, Edwin Barnes & Nicholas J. Mayhall, *Natur. Commun.* 10, 3007 (2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from qiskit_nature.second_q.drivers import PySCFDriver\n",
    "from qiskit_nature.second_q.formats import MoleculeInfo\n",
    "from qiskit_nature.second_q.mappers import ParityMapper\n",
    "from qiskit_nature.second_q.algorithms import GroundStateEigensolver\n",
    "from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD\n",
    "from qiskit_nature.second_q.transformers import FreezeCoreTransformer\n",
    "\n",
    "from qiskit_algorithms.optimizers import SLSQP\n",
    "from qiskit_algorithms.minimum_eigensolvers import AdaptVQE, VQE\n",
    "\n",
    "from qiskit.primitives import Estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Molecule Setup\n",
    "\n",
    "The lithium hydride  at the equilibrium bond length (~1.5 Angstrom) is studied in this tutorial. In order to set up the molecular Hamiltonian one has to get the one- and two-electron integrals from a computational chemistry driver (PySCF in this case) employing the Hartree Fock calculation in STO-3G basis set. The information from the driver is saved in the ``molecule`` object. Then the Fermionic Hamiltonian is mapped to the qubit Hamiltonian using the ``parity`` mapping type as it further reduces the problem size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use PySCF, a classical computational chemistry software\n",
    "# package, to compute the one-body and two-body integrals in\n",
    "# molecular-orbital basis, necessary to form the Fermionic operator\n",
    "\n",
    "info = MoleculeInfo([\"Li\", \"H\"], [(0.0, 0.0, 0.0), (0.0, 0.0, 1.5)])\n",
    "\n",
    "driver = PySCFDriver.from_molecule(info, basis=\"sto3g\")\n",
    "molecule = driver.run()\n",
    "\n",
    "transformer = FreezeCoreTransformer()\n",
    "molecule = transformer.transform(molecule)\n",
    "hamiltonian = molecule.hamiltonian.second_q_op()\n",
    "mapper = ParityMapper(num_particles=molecule.num_particles)\n",
    "tapered_mapper = molecule.get_tapered_mapper(mapper)\n",
    "qubit_op = tapered_mapper.map(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial State and Optimizer\n",
    "Set the initial state to the Hartree Fock state, $| \\psi^{HF} \\rangle$, and choose the optimizer for the VQE algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup the initial state for the variational form\n",
    "init_state = HartreeFock(\n",
    "            molecule.num_spatial_orbitals,\n",
    "            molecule.num_particles,\n",
    "            tapered_mapper,\n",
    "        )\n",
    "\n",
    "estimator = Estimator()\n",
    "\n",
    "optimizer = SLSQP(maxiter=10000, ftol=1e-9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UCCSD-VQE reference\n",
    "Perform the reference calculation employing the ``UCCSD`` variational form and ``VQE`` algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "vqe_ansatz = UCCSD(\n",
    "    molecule.num_spatial_orbitals,\n",
    "    molecule.num_particles,\n",
    "    tapered_mapper,\n",
    "    initial_state=init_state\n",
    ")\n",
    "vqe = VQE(estimator, vqe_ansatz, optimizer)\n",
    "vqe.initial_point = [0] * vqe_ansatz.num_parameters\n",
    "algo = GroundStateEigensolver(tapered_mapper, vqe)\n",
    "result_vqe = algo.solve(molecule)\n",
    "energy_vqe = result_vqe.eigenvalues[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ADAPT-VQE run\n",
    "Now re-instantiate a new ``UCCSD`` object for the use with the ``AdaptVQE`` algorithm (https://qiskit-community.github.io/qiskit-algorithms/stubs/qiskit_algorithms.AdaptVQE.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adapt_vqe_ansatz = UCCSD(\n",
    "    molecule.num_spatial_orbitals,\n",
    "    molecule.num_particles,\n",
    "    tapered_mapper,\n",
    "    initial_state=init_state\n",
    ")\n",
    "\n",
    "adapt_vqe = AdaptVQE(VQE(estimator, adapt_vqe_ansatz, optimizer))\n",
    "result_adapt_vqe = adapt_vqe.compute_minimum_eigenvalue(qubit_op)\n",
    "energy_adapt_vqe = result_adapt_vqe.eigenvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare results\n",
    "Print and compare the gates utilized in both parameterizations and the corresponding energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+---------------------+----------------------------------+\n",
      "|   Ansatz  |   Energy (Hartree)  |              Gates               |\n",
      "+-----------+---------------------+----------------------------------+\n",
      "|   UCCSD   | -1.1001883328961515 | OrderedDict([('EvolvedOps', 1)]) |\n",
      "| ADAPT-VQE | -1.1001795367631413 | OrderedDict([('EvolvedOps', 1)]) |\n",
      "+-----------+---------------------+----------------------------------+\n"
     ]
    }
   ],
   "source": [
    "from prettytable import PrettyTable\n",
    "\n",
    "table = PrettyTable()\n",
    "table.field_names = [\"Ansatz\",\"Energy (Hartree)\",\"Gates\"]\n",
    "\n",
    "vqe_cirq = result_vqe.raw_result.optimal_circuit\n",
    "table.add_row(['UCCSD', str(energy_vqe), vqe_cirq.count_ops()])\n",
    "\n",
    "adapt_vqe_cirq = result_adapt_vqe.optimal_circuit\n",
    "table.add_row(['ADAPT-VQE', str(result_adapt_vqe.eigenvalue), adapt_vqe_cirq.count_ops()])\n",
    "\n",
    "print(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "qiskit==1.1.0\n",
      "qiskit-aer==0.14.2\n",
      "qiskit-algorithms==0.3.0\n",
      "qiskit-ibm-runtime==0.25.0\n",
      "qiskit-machine-learning==0.7.2\n",
      "qiskit-nature==0.7.2\n",
      "qiskit-nature-pyscf==0.4.0\n",
      "qiskit-qasm3-import==0.5.0\n",
      "qiskit-transpiler-service==0.4.5\n"
     ]
    }
   ],
   "source": [
    "! pip freeze | grep qiskit"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
